# include "scene.h"

#ifdef SOLVE_ENERGY
# include <stdio.h>
# include <stdlib.h>

# ifdef _OMP_TEMP_BC_
# include <omp.h>
# endif

extern short interpol_1dim(double const *x, double const *u, unsigned short n, double xp, double *up);
extern double thm_cond(double T);

extern void temp_bc(struct grid *g, double ***T0_imp, short tflag, unsigned short iRK)
   {
    double *r, *z, (*_T)[TEMP_INTRP_N], _z[TEMP_INTRP_N], ***T;
    short i, j, k, nr, nth, nz, _i, _j, tid, nTRD;

   if( (iRK > 0 && tflag != 0) || (iRK == 0 && tflag != 1) || (iRK > 3) )
       {
        fprintf(g->lptr, "\nERROR: In 'temp_bc()': (tflag, iRK) = (%d, %d).\nThis function is designed to handle only following cases:\n\t[tflag = 0, iRK = 1, 2, 3] and [tflag = 1, iRK = 0]\n", tflag, iRK);
        exit(-1);
       }

    nTRD = 1;
#   ifdef _OMP_TEMP_BC_
    nTRD = g->nTRD;
#   endif

    nth = g->nth;

 /* -------------------- Applying BC to the SECTION-0 ------------------------- */
    _T = (double (*)[TEMP_INTRP_N])malloc(nTRD*TEMP_INTRP_N*sizeof(double));

    nr = g->nr[0]; nz = g->nz[0]; r = g->r[0]; z = g->z[0]; T = g->T[0];
    for(_z[0] = z[NzG], k = 1; k < TEMP_INTRP_N; k++) _z[k] = 0.5*(z[NzG+k-1]+z[NzG+k]);

    tid = 0;
#   ifdef _OMP_TEMP_BC_
    #pragma omp parallel private(i,j,k,tid,_i,_j)
#   endif
     {
   /* At the top-boundary open to the upper section. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
         for(k = 0; k < NzG; k++) T[i][j][NzG+nz+k] = g->T[1][i][j][NzG+k];

   /* Dirichlet BC at the INLET. */
#     ifdef _OMP_TEMP_BC_
      tid = omp_get_thread_num();
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+g->nr_nozzle; i++) for(j = NthG; j < nth+NthG; j++)
         {
          _T[tid][0] = T_INLET;
          for(k = 1; k < TEMP_INTRP_N; k++) _T[tid][k] = T[i][j][NzG+k-1];

          interpol_1dim(_z, _T[tid], TEMP_INTRP_N, 0.5*(z[NzG]+z[NzG-1]), &T[i][j][NzG-1]);
         }

   /* At the insulated SOLID BASE after the inlet nozzle. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+g->nr_nozzle; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
         T[i][j][NzG-1] = T[i][j][NzG];

   /* At the insulated cylindrical side-wall. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
        T[NrG+nr][j][k] = T[NrG+nr-1][j][k];

   /* At the top and bottom corner-cells along the cylindrical side-wall. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,2)
#     endif
      for(j = NthG; j < NthG+nth; j++)
          {
           T[NrG+nr][j][NzG-1] = T[NrG+nr-1][j][NzG];
           for(k = 0; k < NzG; k++) T[NrG+nr][j][NzG+nz+k] = g->T[1][NrG+nr][j][NzG+k];
          }

   /* For ghost cells having NEGATIVE radius: applying symmetric B.C. across the center. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,3)
#     endif
      for(k = 0; k < nz+2*NzG; k++) for(j = NthG; j < NthG+nth; j++) for(i = 0; i < NrG; i++)
         T[NrG-i-1][j][k] = T[NrG+i][j][k];

   /* For ghost cells in theta-direction: applying periodic-BC. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,2)
#     endif
      for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
         T[i][NthG-1-j][k] = T[i][NthG+nth-1-j][k];
         T[i][NthG+nth+j][k] = T[i][NthG+j][k];
        }


   /* -------------------- Applying BC to the SECTION-1 ------------------------- */
#     ifdef _OMP_TEMP_BC_
      #pragma omp single
#     endif
       {
// printf("-------\n");
        nr = g->nr[1]; nz = g->nz[1]; r = g->r[1]; z = g->z[1]; T = g->T[1];
       }

   /* The bottom boundary open to the lower section in the central region. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+g->nr[0]; i++) for(j = NthG; j < nth+NthG; j++)
         for(k = 1; k <= NzG; k++) T[i][j][NzG-k] = g->T[0][i][j][NzG+g->nz[0]-k];

   /* At the insulated solid surface between the opening to the bottom section 0 and the extended-outlet section 2. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+g->nr[0]; i < NrG+nr-g->nr_gap; i++) for(j = NthG; j < nth+NthG; j++)
         T[i][j][NzG-1] = T[i][j][NzG];

   /* At the region connected to the extended-outlet 2. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(_i = 0; _i < g->nr_gap; _i++) for(i = _i+NrG+nr-g->nr_gap, j = NthG; j < NthG+nth; j++)
         for(k = 1; k <= NzG; k++) T[i][j][NzG-k] = g->T[2][_i+NrG][j][NzG+g->nz[2]-k];

   /* At the insulated cylindrical side-wall. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
         T[NrG+nr][j][k] = T[NrG+nr-1][j][k];

   /* At the top and bottom corner-cells along the insulated cylindrical solid side-wall */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(j = NthG; j < NthG+nth; j++)
            {
             T[NrG+nr][j][NzG-1] = g->T[2][NrG+g->nr[2]-1][j][NzG+nz-1];
             T[NrG+nr][j][NzG+nz] = T[NrG+nr-1][j][NzG+nz-1];
            }

   /* At the TOP-SURFACE in the central-zone where the heater is embedded. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+g->nr_heater; i++) for(j = NthG; j < nth+NthG; j++)
         T[i][j][NzG+nz] = T[i][j][NzG+nz-1] + (g->t >= t_HEATER_ON ? TOP_HEATER_HF:0)\
                            *(z[NzG+nz]-z[NzG+nz-1])/thm_cond(0.5*(T0_imp[i][j][NzG+nz-1]+T0_imp[i][j][NzG+nz]));

   /* At the insulated solid top-surface "after" the heater-embedded central zone. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,3)
#     endif
      for(i = NrG+g->nr_heater; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
         T[i][j][NzG+nz] = T[i][j][NzG+nz-1];

   /* For ghost cells having NEGATIVE radius: applying symmetric B.C. across the center.*/
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,3)
#     endif
      for(k = 0; k < nz+2*NzG; k++) for(j = NthG; j < NthG+nth; j++) for(i = 0; i < NrG; i++)
         T[NrG-i-1][j][k] = T[NrG+i][j][k];

   /* For ghost cells in theta-direction: applying periodic-BC. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,2)
#     endif
      for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
         T[i][NthG-1-j][k] = T[i][NthG+nth-1-j][k];
         T[i][NthG+nth+j][k] = T[i][NthG+j][k];
        }


   /* -------------------- Applying BC to the SECTION-2 ------------------------- */
#     ifdef _OMP_TEMP_BC_
      #pragma omp single
#     endif
        {
         nr = g->nr[2]; nz = g->nz[2]; r = g->r[2]; z = g->z[2]; T = g->T[2];
        }

   /* The top boundary open to the section-1. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
           T[i][j][NzG+nz] = g->T[1][i+g->nr[1]-g->nr_gap][j][NzG];

   /* At the insulated inner and outer SIDE-WALLs. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
         {
          T[NrG-1][j][k] = T[NrG][j][k]; /* Inner side-wall. */
          T[NrG+nr][j][k] = T[NrG+nr-1][j][k]; /* Outer side-wall. */
         }

   /* Orlansky-BC at the outlet. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,2)
#     endif
      for(_i = 0; _i < nr; _i++) for(i = _i+NrG, j = NthG, _j = 0; _j < nth; _j++, j++)
      switch(iRK)
           {
            case 0:
            T[i][j][NzG-1] = g->outBC.tnext_Tg[_i][_j];
            break;

            case 1:
            T[i][j][NzG-1] = g->outBC.T[2][_i][_j][0]/12 - g->outBC.T[1][_i][_j][0]/2 + 5*g->outBC.T[0][_i][_j][0]/4\
                            + g->outBC.tnext_Tg[_i][_j]/6;
            break;

            case 2:
            T[i][j][NzG-1] = g->outBC.T[2][_i][_j][0]/6 - 3*g->outBC.T[1][_i][_j][0]/4 + g->outBC.T[0][_i][_j][0]\
                            + 7*g->outBC.tnext_Tg[_i][_j]/12;
            break;

            case 3:
            T[i][j][NzG-1] = g->outBC.T[2][_i][_j][0]/24 - g->outBC.T[1][_i][_j][0]/8 + 5*g->outBC.T[0][_i][_j][0]/8\
                            + 11*g->outBC.tnext_Tg[_i][_j]/24;
            break;
           }

   /* For the corner-cells. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,1)
#     endif
      for(j = NthG; j < NthG+nth; j++)
          {
           T[NrG+nr][j][NzG+nz] = g->T[1][NrG+g->nr[1]-1][j][NzG];
           T[NrG-1][j][NzG+nz] = g->T[1][NrG+g->nr[1]-g->nr_gap-1][j][NzG];

        /* These require application of the outlet-BC first. */
           T[NrG-1][j][NzG-1] = T[NrG][j][NzG-1];
           T[NrG+nr][j][NzG-1] = T[NrG+nr-1][j][NzG-1];
          }

   /* For ghost cells in theta-direction: applying periodic-BC. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(dynamic,1)
#     endif
      for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
         T[i][NthG-1-j][k] = T[i][NthG+nth-1-j][k];
         T[i][NthG+nth+j][k] = T[i][NthG+j][k];
        }

     }

    free(_T);

    return;
   }
#endif
